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ABSTRACT 

In this paper we discuss the influence of gravitational instabilities in massive pro¬ 
tostellar discs on the dynamics of dust grains. Starting from a Smoothed Particle 
Hydrodynamics (SPH) simulation, we have computed the evolution of the dust in a 
quasi-static gas density structure typical of self-gravitating disc. For different grain size 
distributions we have investigated the capability of spiral arms to trap particles. We 
have run 3D radiative transfer simulations in order to construct maps of the expected 
emission at (sub-)millimetre and near-infrared wavelengths. Finally, we have simu¬ 
lated realistic observations of our disc models at (sub-)millimetre and near-infrared 
wavelengths as they may appear with the Atacama Large Millimetre/sub-millimetre 
Array (ALMA) and the High-Contrast Coronographic Imager for Adaptive Optics (Hi- 
CIAO) in order to investigate whether there are observational signatures of the spiral 
structure. We And that the pressure inhomogeites induced by gravitational instabil¬ 
ities produce a non-negligible dynamical effect on centimetre sized particles leading 
to significant overdensities in spiral arms. We also And that the spiral structure is 
readily detectable by ALMA over a wide range of (sub-)millimetre wavelengths and 
by HiCIAO in near-infrared scattered light for non-face-on discs located in the Ophi- 
ucus star-forming region. In addition, we And clear spatial spectral index variations 
across the disc, revealing that the dust trapping produces a migration of large grains 
that can be potentially investigated through multi-wavelenghts observations in the 
(sub-)millimetric. Therefore, the spiral arms observed to date in protoplanet ary disc 
might be interpreted as density waves induced by the development of gravitational 
instabilities. 

Key words: accretion, accretion discs - gravitation - instabilities - stars: pre-main- 
sequence - sub-millimetre: stars - planet and satellites: formation. 


1 INTRODUCTION 

The physical and chemical evolution of protostellar 
discs plays a crucial role in the planet formation process. 
According to the most widely accepted scenario of planet for¬ 
mation, the “core accretion” model, in the high-density envi¬ 
ronment of a circumstellar disc, micron sized dust grains are 
expected to grow via collisional agglomeration to kilometre¬ 
sized planetesimals (Testi et al. 2014). In recent years, the 
investigation on the dust evolution in protostellar discs has 
been the subject of several theoretical (e.g. Brauer et al. 
2008, Birnstiel et al. 2010) and laboratory studies (e.g. Blum 
et al. 2000, Blum and Wurm 2008) aimed at understand- 
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ing dust dynamics and its growth. The general picture that 
emerges from these investigations is that the growth from 
submicron sized particles to larger objects is a complex 
process that contains many physical challenges (Zsom and 
Dullemond 2008, Okuzumi 2009, Birnstiel et al. 2010, Laibe 
et al. 2012). 

Since the dust-to-gas mass ratio in protostellar discs is 
much lower than unity (typical value of ~ 10“^), the dynam¬ 
ics of dust particles is heavily affected by the aerodynamic 
drag force that arises from the velocity difference between 
dust particles and the surrounding gas. For standard disc 
geometries, the pressure in the circumstellar disc tends to 
decrease outwards leading the gas to orbit at sub-Keplerian 
velocities. On the other hand, since the dust is not affected 
by the gas pressure, a freely orbiting dust particle feels only 


© 0000 RAS 



2 Giovanni Dipierro, Paola Pinilla, Giuseppe Lodato and Leonardo Testi 


centrifugal forces and gravity, and therefore, orbits at Ke- 
plerian velocity to a first approximation. As a result, the 
drag force removes the dust angular momentum leading to 
an inward drift at a rate that depends on the size of par¬ 
ticles (Weidenschilling 1977). The radial velocity induced 
by gas drag is high enough to produce a rapid depletion of 
mm-sized particles in the outer disc before they can become 
large enough to decouple from the gas and, thus, preventing 
the growth of planetesimals required for the formation of the 
planetary cores. This is what is known as radial drift barrier, 
and it still represent one of the main unsolved problem of 
the early phases of planet formation (Weidenschilling 1977, 
Nakagawa et al. 1986, Klahr and Henning 1997, Brauer et al. 
2008). 

In addition, the growth of solid particles has to circum¬ 
vent another hurdle closely linked to the sticking efficiency of 
colliding dust grains. It is expected that, while low-velocity 
collisions result in grain growth, high-velocity impacts be¬ 
tween particles are destructive (Poppe et al. 1999, Blum 
et al. 2000, Blum 2006). From laboratory and numerical 
studies it has been shown that, as particles become larger, 
they tend to collide at higher impact velocities. Therefore, 
it is expected that the velocities of mm-sized objects in the 
outer disc are larger than the critical threshold for stick¬ 
ing, so that collisions lead to fragmentation and thus pre¬ 
vent dust particles from forming larger bodies (Brauer et al. 
2008). As a result, the growth of dust particles toward meter 
sizes needs to overcome two obstacles: the rapid depletion 
of material due to the radial drift and the particle fragmen¬ 
tation produced by destructive collision. 

From an observational point of view, the evidence that 
grains in discs are significantly larger than those in the dif¬ 
fuse interstellar medium (ISM) has been obtained from a va¬ 
riety of observational techniques (see the reviews by Natta 
et al. 2007 and Testi et al. 2014). Multi-wavelenght (sub-) 
mm observations (e.g., Testi et al. 2001, 2003, Ricci et al. 
2010, Guilloteau et al. 2011, Perez et al. 2012) of protostel- 
lar discs have shown that grain growth is efficient enough 
to quickly produce mm-sized particles which are retained in 
the outer disc for a relatively long time. Therefore, how to 
theoretically explain the retention of mm-sized particles in 
the outer regions of disc remains an open question. 

Several ideas of plausible mechanisms aimed at slowing 
down the rapid inward drift have been proposed, such as vor¬ 
tices, planet-disc interactions, snowlines, zonal flows or dead 
zones (e.g. Klahr and Henning 1997, Youdin and Shu 2002, 
Johansen et al. 2009, 2011, Pinilla et al. 2012). The common 
feature of these mechanisms is to create long-lived radial 
and/or azimuthal inhomogeneites in the gas density struc¬ 
ture in order to “trap” particles into the overpressure regions 
as discussed initially by Haghighipour and Boss (2003a,6). 
Some observational support for particle trap has been ob¬ 
tained from recent high angular resolution observations at 
mm-wavelength using ALMA that suggested the presence of 
a particle trap by an anti-cyclonic vortex that creates dust 
segregation of small and large particles, as for example the 
case of IRS 48 (van der Marel et al. 2013). 

In this paper we discuss the influence of the spiral den¬ 
sity waves induced by gravitational instabilities (GI) on the 
dust dynamics. It is believed that in the early stage of star 
formation (Class 0/Class I objects), the disc could be mas¬ 
sive enough to have a non-negligible dynamical effect on the 


evolution of the overall system. The disc self-gravity may 
affect the disc dynamics through the propagation of den¬ 
sity waves that lead to the formation of a prominent spiral 
structure (Lodato and Rice 2004, 2005, Durisen et al. 2007, 
Cossins et al. 2009, Forgan et al. 2011). In this context, it 
is expected that the density inhomogeneities induced by the 
development of GI affect both the dynamics and the growth 
of the dust grains. Recent work has already suggested that 
grain growth may be accelerated by concentrating the grains 
in the peak of the self-gravitating structures, possibly over¬ 
coming the issue of the rapid grain migration (Rice et al. 
2004, 2006, Gibbons et al. 2012, 2014). Moreover, the dust 
retention in spiral arms implies low relative velocities be¬ 
tween particles, which suggests that collisions are more likely 
to be constructive (Gibbons et al. 2012, 2014). 

In addition, the dust migration towards spiral arms 
can produce macroscopic features that can be identified 
through high angular resolution observations. Spiral and 
non-axisymmetric structures in protostellar discs have been 
observed recently in transitional discs by imaging the distri¬ 
bution of scattered light at near-infrared wavelengths (Muto 
et al. 2012, Garufi et al. 2013, Grady et al. 2013) and by 
high-resolution ALMA observations of molecular line emis¬ 
sion (Fukagawa et al. 2006, Ghristiaens et al. 2014). For 
such evolved discs, one generally expects the disc not to be 
massive enough to be self-gravitating, and thus the more 
common explanation for the origin of the spiral is the dy¬ 
namical interaction with an embedded planet. Nonetheless, 
it is worth remarking that disc mass estimates suffer from 
systematic errors, mostly due to uncertainties in the dust 
opacity and the assumption of a constant dust-to-gas ratio. 
Therefore, it is still debatable if these observed spirals are 
or not produced by planet-disc interactions (Juhasz et al. 
2014). Moreover, Dipierro et al. (2014) (see also Gossins 
et al. 2010) have found that, assuming that the dust is per¬ 
fectly mixed with the gas, spiral structures of a variety of 
non-face-on self-gravitating circumstellar discs models with 
different properties in mass and radial extension are read¬ 
ily detectable by ALMA over a wide range of wavelengths. 
One of the aims of the present work is to improve the ob¬ 
servational predictions of self-gravitating circumstellar discs 
presented in Cossins et al. (2010) and Dipierro et al. (2014) 
by including the treatment of dust dynamics (in particular, 
by relaxing the hypothesis that gas and dust are perfectly 
mixed) and a more detailed procedure for the generation of 
emission maps. 

In this work, by combining the results of hydrodynamic 
dust and gas simulations with 3D Monte Carlo radiative 
transfer calculations (using RADMC-3D), we study if the 
occurrence of local pressure maxima induced by GI can trap 
dust particles and investigate the detectability of these in¬ 
homogeneites at near-infrared and (sub-)milfimetre wave¬ 
lengths. 

The paper is organized as follows: in Section 2 we de¬ 
scribe the details of the hydrodynamic simulations. In Sec¬ 
tion 3 we describe the results of the simulations of the dust 
dynamics and present the observational predictions of our 
disc models. Finally, in Section 4 we discuss the significance 
of our results. 
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Figure 1. Gas surface density structure of a O.25M0 disc around 
a IMq star that has reached a long-lived, self-gravitating state. 



Figure 2. Gas volume density, temperature and pressure gradient 
along a slice across the disc. 


2 GAS AND DUST DISC MODELS 

The simulations presented here calculate the evolution 
of the dust surface density in a quasi-static gaseous self- 
gravitating disc. The system comprises two components 
discs: a gas disc, characterized by a spiral surface density 
structure induced by GI, and a dusty disc that evolves un¬ 
der the action of drag forces and turbulent mixing induced 
by the gas-dust aerodynamical coupling. The self-gravitating 
gas disc model has been taken from a snapshot of a three- 
dimensional Smoothed Particle Hydrodynamics (SPH) sim¬ 
ulation performed by Lodato and Rice (2004). This model 
can be considered as representative of a self-gravitating gas 
disc in a quasi-steady state characterized by a spiral sur¬ 
face density structure produced by a combined effect of self¬ 
gravity and cooling. We compute the dust surface density 
evolution in such fixed gaseous background in the radial di¬ 
rection along different azimuthal cuts across the disc. After¬ 
wards, starting from the resulting surface density structure 
we create 3D models using a representative estimator for the 
tickness of the dusty disc. We carry out 3D radiative transfer 
simulations in order to create the expected emission maps at 
(sub-)millimetre and near-infrared wavelengths of our mod¬ 
els which are used to simulate realistic observations. 

2.1 Gas disc model 

The gas disc model used as an input for the dust simu¬ 
lations has been taken from Lodato and Rice (2004). They 
have performed three-dimensional SPH simulations of self- 
gravitating discs evolving under an analytically prescribed 
cooling function. The system comprises a central star mod¬ 
elled as a point mass M* surrounded by a gaseous disc of 
mass Mdisc which extends from Rin — 0.25 to Rout = 25 
in code units. The simulation used in this work was per¬ 
formed using a value for the disc-to-central object mass ratio 
q = Mdisc/Af* = 0.25. The disc is characterized initially by 
a surface density profile S oc R~^ while the initial tempera¬ 
ture profile is T oc Regarding the thermal aspect of 

disc evolution, the heating is governed by both PdV work 
and viscous dissipation while the cooling is implemented us¬ 


ing a simple cooling law given by (Gammie 2001): 


where Ui is the specific internal energy and tcooi.i is the cool¬ 
ing time associated with the particle. The latter is de¬ 
termined using the simple parametrization Ditcooi.i = /3cooi, 
where Pcooi is held constant and equal to 7.5. All simula¬ 
tions described above are essentially scale-free. In order to 
consider a gas disc model with realistic properties in mass 
and radial extension, we modify the length and the mass 
scale of the data output of the SPH simulations such that 
the radial extent of the disc is 150 au with a disc mass equal 
to 0.25 Mq. As shown in Lodato and Rice (2004), the ini¬ 
tial phase of disc evolution is governed by a rapid cooling 
until the gravitational instability becomes effective, provid¬ 
ing a source of effective heating through both compression 
and shocks. In approximately one thermal time, the disc 
settles into a quasi-steady state developing spiral density 
waves that propagate across the disc. Once the thermal equi¬ 
librium state is reached, the gas disc is characterized by a 
spiral structure (shown in Fig. 1) where the heating gener¬ 
ated through GI balances the external cooling. This struc¬ 
ture remains quasi-steady, evolving secularly on the viscous 
timescale (10®-10® yrs in this case). Individual spiral arms 
are recurrent structures evolving on the dynamical timescale 
(~ 10® yrs in the outer disc). 

Fig. 2 shows the volume density, temperature and pres¬ 
sure gradient along a slice across the disc. It can be noticed 
that all the gradients are characterized by a similar trend 
and are null at the same radial positions. Therefore, the 
overdensity regions, which correspond to the arm regions, 
are characterized by relatively higher temperatures while in 
the interarm region the temperature is lower (see also Fig. 
3 of Dipierro et al. 2014). This is due to the increase in 
the pressure caused by a shock driving a compression of gas 
and, as a consequence, a growth of the temperature. In the 
interarm regions, the increase in the pressure provides an ex¬ 
pansion of gas parallel to the shock and, as a consequence, a 
decrease of the temperature (Gossins et al. 2009). Moreover, 
since the pressure depends on the volume density and the 
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sound speed, the peaks of the pressure gradient are higher 
than those related to the temperature and volume density 
profile. 


2.2 Dust evolution model 

To describe the evolution of the dust due to drag forces 
and turbulent diffusion it is convenient to define how well 
coupled the particles are to the gas. The gas-dust aerody¬ 
namical coupling is expressed by the Stokes number, usually 
denoted by St, defined as the ratio of the stopping time of a 
particle and the local dynamical timescale Generally, 

there are four different regimes for the dust-to-gas aerody¬ 
namical coupling depending on the value of the Reynolds 
number and the ratio between mean free path Amfp of the 
gas molecules and the dust particle size, a. In the Epstein 
regime, for Amfp/a 1? 4/9, the Stokes number for particles 
close to the disc mid-plane, is given by (Birnstiel et al. 2010): 


Generally, it is expected that a gravitationally unsta¬ 
ble disc cannot be treated as a viscous disc because gravity 
can affect the transport mechanism over long distances. In 
this context, a number of numerical simulations have been 
recently carried out with the aim to test the locality of trans¬ 
port provided by GI. Lodato and Rice (2004), Cossins et al. 
(2009) and Forgan et al. (2011), using an a-like approach, 
have compared the dissipation power provided by gravita¬ 
tional instability with the expectations based on a viscous 
theory of disc. They have found that for disc to star mass 
ratio Mdisc/M* 0.25, the gravitational instability is well 
described using a viscous approach. 

Encompassing all the systematic and random motions 
described above, the time evolution for the dust surface den¬ 
sity EJj of the i*** dust species with size a‘ can be described 
by an advection-diffusion equation (Brauer et al. 2008, Birn¬ 
stiel et al. 2010): 


ask 1 

dt RdR 




= 0 , 


( 5 ) 


St = 


aps n 
^2 ’ 


( 2 ) 


where Eg is the gas surface density and ps is the internal 
density of dust grains. The radial velocity of dust particles 
can be determined by self-consistently solving the radial and 
azimuthal components of the momentum equation (Weiden- 
schilling 1977): 


Ur — Udrag “t“ Udrift 


Mg,r 1 1 dP 

1 -I- St^ St -I- St“^ Pgfl dR ’ 


( 3 ) 


where fl is the Keplerian angular velocity, P is the pressure 
and pg is the gas volume density. The first term is the gas 
drag term which is produced by the coupling of the dust 
particle to gas moving with a radial velocity of Mg,r given 
by viscous evolution. The second term is the radial drift 
caused by the angular momentum exchange between dust 
particles and the non-Keplerian orbital motion of the gas. If 
the Stokes number is much smaller than unity, the dust par¬ 
ticle is strongly coupled to the gas, which implies that the 
gas and dust share the same motion. On the other hand, 
for St >> 1, the particles are unaffected by the gas and 
consequently move on Keplerian orbits to a first approxi¬ 
mation. Particles in the intermediate size range, i.e. St ~ 1, 
have large drift velocities and are expected to experience the 
highest concentration in pressure maxima. 

Since the gas is expected to be turbulent, the dust is 
mixed by the gas due to aerodynamical coupling. The tur¬ 
bulent transport depends on the gas diffusion Dg, which 
is usually described by the canonical turbulence prescrip¬ 
tion (Shakura and Sunyaev 1973). The dust diffusivity in¬ 
duced by the gas-dust aerodynamical coupling is defined as 
(Youdin and Lithwick 2007): 


D 
1-f St^ 


a{cs)Hg 

l-fSt^ ’ 


( 4 ) 


where, as previously described, the diffusion coefficient (see 
Eq. 4) and the radial velocity (see Eq. 3) depends on the 
Stokes number and therefore on the particle size. 

In this work, the dust evolution for each grain size in our 
sample is computed using the numerical method adopted in 
Birnstiel et al. (2010). All the simulations performed here are 
evolved for one outer orbital period: 10^ yrs. Such timescale 
is long enough to ensure that the drag forces and turbulent 
diffusion affect the dust density structure but short enough 
that we do not expect a significant evolution of the spiral 
structure. We cut the disc into 720 radial directions. For 
each cut, we compute the gas pressure needed in Eq. 3. We 
then solve, for each cut, Eq. 5 using 100 radial cells, so as 
to obtain the 2D surface density of the dust. 


2.3 Initial conditions 

The evolution of the dust density and the resulting ob¬ 
servational signatures depend on the initial conditions of 
the simulations. Generally, the distribution of grain size in 
a protostellar disc is affected by the grain growth and frag¬ 
mentation processes which in turn are closely linked to the 
relative velocities between colliding particles. In this work, 
we assume that the dusty disc has settled into a coagulation- 
fragmentation steady state. Since the self-gravitating phase 
of protostellar disc is more likely to occur in the early stage 
of star formation, such initial condition for the dust density 
ditribution might appear simplistic. However, recent obser¬ 
vational results have shown that large dust aggregates can 
form during the disc formation stage in the infalling enve¬ 
lope (Miotello et al. 2014). The distribution of grain size is 
assumed to be a standard power-law size distribution: 

n(a) OC a~‘‘ for Omin < a < Umax , (6) 


where Hg is the gas scale height, (cs) is the azimuthal aver¬ 
age of the sound speed and a is an unknown dimensionless 
parameter which embodies all the uncertainties about the 
nature of the viscous mechanism. Since the diffusion time 
scale is longer than the lifetime of individual spiral features 
in self-gravitating discs (~ D~^), the dust diffusivity is com¬ 
puted using the azimuthal average of the sound speed. 


where Umin and Umax represent the minimum and maximum 
size of the grains. The dust size distribution is modelled 
using 100 sizes and it is normalized such that the total initial 
dust density is one per cent of the gas density in each cell of 
the grid. The minimum size of the grains is assumed to be 
equal to 0.1 pm. 

Regarding the value of Umax, Birnstiel et al. (2010) have 
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found that an acceptable estimator for the maximum grain 
size of the dust distributions can be computed by assuming 
the random turbulence motion as the only source of relative 
motions (for a detailed discussion of the turbulent relative 
velocities see Ormel and Cuzzi 2007). In this context, Rice 
et al. (2005) have found that the turbulent diffusion induced 
by GI for the gas disc model used here, is characterized by 
a “gravoturbolent” parameter a = 0.05. Since the turbu¬ 
lence parameter related to GI is reasonably larger than the 
typical turbulence parameters, the previously mentioned as¬ 
sumption for the maximum grain size of the dust size distri¬ 
bution can be considered a valid approximation. The value 
of the maximum grain size is therefore computed by equat¬ 
ing the relative velocities between colliding particles induced 
by “gravoturbolent” diffusion and the threshold velocity for 
fragmentation (Birnstiel et al. 2010). In the Epstein regime, 
this is given by: 


_ 4(Sg) 't^frag 
STvaps {cs)'^ 


(7) 


where wnag is the fragmentation threshold velocity and (Eg) 
is the azimuthal average of the gas surface density. Similarly 
to the computation of the dust diffusivity (Eq. 4), since the 
typical grain growth timescales is longer than the lifetime 
of individual spiral features in self-gravitating discs (Brauer 
et al. 2008, Okuzumi et al. 2011), Umax is computed using 
the azimuthal average of the gas surface density and sound 
speed. Regarding the value of the fragmentation velocity, it 
has been shown through theoretical and experimental stud¬ 
ies that it depends on dust grains material properties. Typi¬ 
cally, the fragmentation velocity for silicates dust aggregates 
is of the order of a few m s“^ (Blum and Wurm 2008) while 
icy aggregates are able to grow at collisions with velocities up 
to 50 m s“^ (Wada et al. 2009, Gundlach and Blum 2015). 

Moreover, the dust size distribution is characterized by 
the value of the power-law exponent q that depends on the 
dust dynamics and the processes of grain coagulation and 
fragmentation that occur in the disc. While the typical value 
in the ISM is qism = 3.5 (Mathis et al. 1977), in a protoplan¬ 
etary disc if coagulation (fragmentation) processes dominate 
over fragmentation (coagulation), a lower (higher) value for 
q than the ISM one is expected. 

In order to cover a wide range for the initial dust size 
distributions, we perform simulations using different values 
for the fragmentation threshold velocity and the power-law 
exponent. We consider two values for the fragmentation 
threshold velocity: unag = [10,30]ms“^ and three values 
for the power-law exponent: q = [3.0, 3.5,4.0]. 


2.4 Monte Carlo radiative transfer simulation 

We calculate the expected emission maps of the self- 
gravitating disc model adopted here using the 3D radiative 
transfer code RADMG-3D^. We focus on images at H-band 
(1.65 /rm) and ALMA band 6 (220 GHz), band 7 (345 GHz), 
band 8 (460 GHz) and band 9 (680 GHz). Our aim is to 
compute HiGIAO scattered light images in H-band polar¬ 
ized intensity and ALMA (sub-)millimetre images in order 


^ http://www.ita.uni-heidelberg.de/ dullemond/software/radmc- 
3d 


to study the detectability of the peculiar spiral structure in¬ 
duced by GI both in scattering and thermal dust emission. 

It is known that the scattering intensity is closely linked 
to the behaviour of the scattering phase function of dust 
grains. While small grains scatter photons uniformly, the 
scattering phase function of large grains is dominated by for¬ 
ward scattering. Therefore, since large grains scatter stellar 
photons mostly into the disc, tracing larger grains by analyz¬ 
ing scattering images of non-edge-on discs results challeng¬ 
ing (Mulders et al. 2013). Thus, polarized scattering images 
at short wavelengths trace small (0.1-10 pm) dust grains at 
the surface of the disc (where stellar photons get absorbed or 
scattered). At (sub-)millimetre wavelengths, since the disc 
is generally optically thin in the vertical direction, the emis¬ 
sion probes the mid-plane disc surface density in large grains 
(0.1-10 mm), due to their high opacity at these wavelengths 
(Dullemond et al. 2007, Williams and Cieza 2011). 

In order to perform RADMG-3D simulation, we need 
to create 3D structures starting from the 2D dust density 
distributions obtained with the procedure described in Sect. 
2. In a protostellar disc the vertical structure of the dust is 
determined by the balance between the dust settling toward 
the disc mid-plane due to the gravity and the vertical tur¬ 
bulent stirring induced by the gas drag. The dust density 
profile of the i‘*' dust species with size a' is assumed to be 
Gaussian: 


Pd(z) = 




( 8 ) 


where the scale height is computed assuming that 
the gravitational settling balances the turbulent diffusion 
(Brauer et al. 2008, Birnstiel et al. 2010): 


h‘^ — HglR) min 1, 


min(St‘, 1/2)(1 + St' ) 


( 9 ) 


where, bearing in mind that the disc is in margina self- 
gravitating state, the gas scale height is given by (Gossins 
et al. 2009): 


Bg(R) 


7TEg(R)R^ 

M* 


( 10 ) 


Having constructed the 3D distribution for all the dust 
species in our sample, we compute the temperature struc¬ 
ture of the dust via thermal Monte Garlo simulations. The 
source of radiation is assumed to be the central star, located 
at the centre of the coordinate system, with M* = Mq, 
Ri, = Rq and TeS = 6000 K. The dust in our model consists 
of spherical composite porous grains with chemical abun¬ 
dances adopted in Ricci et al. (2010). The absorption and 
scattering cross sections for each dust species in our sample 
are calculated using the routine developed by Bohren and 
Huffman (1983). In addition, in order to compute near in¬ 
frared scattered light images in the H-band, we include the 
full treatment of polarization by computing the scattering 
matrix for each dust species. 

Starting from the dust temperature structure of the 
model, the expected emission maps are computed via ray¬ 
tracing assuming that the disc is non-face-on (inclination an¬ 
gle of 30°). The polarized intensity image (PI=\/Q^”+TP) 
is calculated starting from the full Stokes vector (I,Q,U,V) 
computed by RADMC-3D. The full-resolution (sub-)mm im- 
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ages directly produced by RADMC-3D are used as input sky 
models to simulate realistic ALMA observations using the 
Common Astronomy Software Application (CASA) ALMA 
simulator. The HiCIAO images are computed by convolving 
the full-resolution images in H-band with a measured Hi¬ 
CIAO point spread function taken from the ACORNS-ADI 
SEEDS Data Reduction Pipeline software (Brandt et al. 
2013). 


3 RESULTS 

In this section, we present results of the dust simulations 
and the observational predictions of the resulting disc mod¬ 
els. We will focus on the level of the arm-interarm contrast 
induced by gas-dust coupling for small and large grains and 
how the dust evolution and observational signatures change 
with different initial grain size distributions. 


3.1 Dust density structure 

First, we show a representative set of results with a 
given set of model parameters. Figure 3 shows the initial 
(left) and the final (right) dust density distribution along a 
slice across the disc from a simulation starting from an initial 
size distribution with exponent q = 3.5 and a fragmentation 
velocity equal to 30 m s“^. The colour contours denote the 
dust surface density distribution o^a, R) per logarithmic size 
bin, defined through: 

/ ®max (.R) 

cr{a, R) dlna , (11) 

min 

where Ed is the total dust surface density. The solid line 
shows the fragmentation barrier, i.e. the maximum value of 
the grain size distribution (see Eq. 7). The dashed line rep¬ 
resents the particle size corresponding to a Stokes number 
of unity, which represent the size of particles that experi¬ 
ence the fastest radial drift. Since the dust is in the Epstein 
regime, this line reflects the shape of the surface density, 
as expressed in Eq. 2 so that peaks in this quantity refer 
to the high density arms, while valleys indicate the inter¬ 
arm regions. The dotted line in Fig. 3 shows the particle 
size corresponding to a Stokes number of a^Pcooi, which, 
as described below, represent the size above which the ra¬ 
dial migration towards pressure maximum dominates over 
turbulent diffusion. It can be noticed that particles with 
Stokes number St ~ 1 tend to be affected by spiral density 
waves most effectively and exhibit the largest concentration 
in wave crests. The dust migration toward pressure maxima 
is negligible for smaller particles. This is due to the higher 
coupling between the dust and the slowly evolving gas and 
also because the diffusion is strong enough such that vari¬ 
ations in the dust-to-gas ratio are being smeared out. In 
addition, since the maximum value of the initial grain size 
distribution decreases with radius, there is a strong den¬ 
sity gradient of the largest particles between adjacent cells 
along the radial grid. Therefore, the gravoturbulent diffu¬ 
sion produces a fast migration of larger particles towards 
outer regions leading to an additional density enhancement 
in the spiral arm as can be easily noticed in the right panel 
of Fig. 3, where a large population of particles above the 
fragmentation limits is present. The motion towards outer 


regions induced by diffusion is more evident when the slope 
of maximum grain size increases (for R > 70 au) because 
a higher number of particles diffuse towards the outer ra¬ 
dius. However, it is expected that largest particles beyond 
the fragmentation threshold tend to experience disruptive 
collisions on longer timescales. 

It is known that the dust concentration in pressure 
bumps is the result of the competition between the radial 
drift, that forces the particle to migrate towards pressure 
maxima, and the turbulent diffusion that tends to recover 
the dust density distribution into the initial condition where 
the dust-to-gas ratio is equal in the overall disc. In order to 
hgure out the relative importance of advection and diffusion, 
we can compare the advection time scale, that can be esti¬ 
mated by the minimum time required to concentrate solid 
material in a pressure maximum: tadv ~ A/udrift, and the 
diffusion time scale tdifi ~ A^/Dd, where A is the scalelength 
of the density inhomogeneities that, in a gravitationally un¬ 
stable disc, is of the order of the disc scale height H^. The 
ratio of these timescale is given by: 

^diff _ St ASg , . 

, ... 

tadv 

where, since we have a surface density enhancement ~ 
ASg over a radial length scale Hg, we approximate 
{1/pg)(dP/dR) as ~ (Cs/I7g)(AEg/Sg). Thus, for particles 
with St L a (ASg/Eg)“^, the radial migration towards pres¬ 
sure maximum dominates over turbulent diffusion. In a self- 
gravitating disc in thermal equilibrium, the amplitude of the 
spiral structure induced by gravitational instability has to 
provide a stress large enough to balance the effective cool¬ 
ing rate (see Eq. 1). Thus, it is expected that the strength of 
the modes are proportional to the cooling rate. In this con¬ 
text, by performing high-resolution three-dimensional SPH 
simulations of self-gravitating discs evolving under the an¬ 
alytically prescribed cooling function expressed by Eq. 1, 
Cossins et al. (2009) have found that the strength of the 
modes varies such that: 


AEg 


VIP' 


(13) 


where (AEg/Eg) indicates the relative RMS amplitude of 
the surface density perturbations. This is closely linked to 
the fact that smaller cooling times result in larger pres¬ 
sure gradients and thus a more efficient dust trapping into 
the pressure maxima (Gibbons et al. 2012). Moreover, the 
value of the gravoturbulent parameter a depends on the 
imposed cooling. In detail, using the general results of vis¬ 
cous disc theory (see Section 2 of Lodato and Rice 2004), 
the value of a needed to balance the imposed cooling varies 
such that a oc Thus, the minimum Stokes number 

above which the advection towards pressure maxima dom¬ 
inates over gravoturbulent diffusion varies with the cooling 
as . In other words, the range extension of particle size 

that show the concentration into pressure maxima decreases 
with decreasing cooling time due to the combined effect of 
radial drift and turbulent mixing. 

In Fig. 4 we show the final surface density rendered im¬ 
ages of the dust for all whole disc. In order to show the dif¬ 
ferential distribution from different ranges of dust grain size, 
the upper panels show the surface density integrated in size 
in the range [1/rm, 100/rm] (left) and [1cm, lOOcm] (right). 
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R [au] R [au] 

Figure 3. Dust density distribution at the beginning (left) and at the end (right) of simulation {t = 10^ yrs) along a slice across the 
disc. The simulation starts from an initial size distribution with exponent q = 3.5 and a fragmentation velocity equal to 30 m s“^. The 
dashed line represents the particle size corresponding to a Stokes number of unity while the dotted line show the particle size where 
St = (yy/0coo\- The solid line shows the fragmentation barrier, i.e. the maximum value of the initial grain size distribution (see Eq. 7). 
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Figure 4. The upper panels show the surface density at the end of simulation {t = 10^ yrs) integrated in size in the range [Ifim, lOOfim] 
(left) and [1cm, 100cm] (right). The lower panels show the total dust density distribution (left) and the dust-to-gas ratio (right) at the 
end of simulation. 
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Figure 5. Dust density distribution at the beginning (left) and at the end (right) of simulation (t = 10^ yrs) along a slice across the 
disc. The simulation starts from an initial size distribution with exponent q = 3.5 and a fragmentation velocity equal to 10 m s“^. The 
dashed line represents the particle size corresponding to a Stokes number of unity while the dotted line show the particle size where 
St = 0\//5cool- The solid line shows the fragmentation barrier, i.e. the maximum value of the initial grain size distribution (see Eq. 7). 


As expected, the dust density distribution for particles with 
Stokes number close to unity (see the upper right panel of 
Fig. 4) show the highest concentration in spiral arms. In this 
case, the spiral arms are very thin and display a strong con¬ 
trast between arm and Interarm regions, while the structure 
shown by smaller particles is broader and essentially follows 
the gas structure (see Fig. 1). Note that beyond ~ 50 au the 
density of large particles drops because we have crossed the 
fragmentation threshold. The lower panels of Fig. 4 show the 
total dust density distribution (left) and the dust-to-gas ra¬ 
tio (right). It can be noticed that the dust transport process 
in a self-gravitating accretion discs quickly modify the radial 
and azimuthal distribution of the dust-to-gas ratio. Starting 
from an initial uniform dust-to-gas ratio (10“^), dust parti¬ 
cles migrate quickly towards pressure maxima leading to a 
creation of a spiral structure in the dust-to-gas ratio. The 
local density of particles in interarm regions reaches levels 
a factor of ~ 10~^ times that of the initial dust density. It 
is worth noting that the evolution of the dust-to-gas ratio 
shows a remarkable depletion of dust grains from the inter¬ 
arm regions since larger particles are more subject to radial 
drift due to the lower gas density in interarm regions that 
produce a better coupling between dust and gas (see Eq. 2). 

3.1.1 Influence of fragmentation velocity 

In order to cover a wide range of initial dust density dis¬ 
tribution, we run the same simulations, but we now change 
the value of the fragmentation threshold velocity. As men¬ 
tioned above, dust aggregates are able to grow by collisions 
with velocities up to a value that depends on the material 
properties (Blum and Wurm 2008, Wada et al. 2009). 

In Fig. 5 we show the initial (left) and the final (right) 
dust density distribution along a slice across the disc using a 
fragmentation velocity equal to 10 m s“^ and an initial size 
distribution with exponent q = 3.5. If the velocity at which 
particles fragment decreases from 30 to 10 m s“^, grains are 
not expected to grow to large sizes so that the condition 



Figure 6. Dust-to-gas mass ratio at the end of the simulation 
{t = 10^ yrs) along a slice across the disc for two values of the 
power-law exponent of the grain size distribution: q = 3,4. The 
dotted line shows the initial dust-to-gas ratio. 

St ~ 1 is not expected to be met anywhere in the disc. In 
this case, the efficiency of the dust trapping is reduced since 
most of the particles are characterised by a Stokes number 
less than unity. As a result, variations in the dust-to-gas 
ratio are being smeared out by the turbulent diffusion and 
the radial drift has only a minor effect in the dust dynamics. 

3.1.2 Influence of the initial power-law grain size 
distribution 

As mentioned above, the grain coagulation and frag¬ 
mentation processes can affect the size distribution quite 
significantly. If the fragmentation process dominates over co¬ 
agulation, a constant replenishment of smaller grains occurs 
in the disc producing a steeper size distribution, approaching 
n{a) oc a“^. In this case, most of the mass will be concen- 
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trated in small grains which are less affected by the radial 
drift toward pressure maxima. By contrast, if the grain co¬ 
agulation dominates over fragmentation, the resulting grain 
size distribution is expected to be Hatter (e.g. n(a) oc a~^). 
In this case, the grain coagulation process would produce 
large grains which migrate quickly toward pressure maxima. 
Fig. 6 shows the radial dependence of the final dust-to-gas 
mass ratio along a slice across the disc for two values of the 
power-law exponent of the grain size distribution and using 
a value of the fragmentation velocity equal to 30 m s“^. For 
g = 3, the radial distribution of dust is strongly affected by 
the migration toward pressure maxima induced by gas drag 
since most of the dust mass is in large grains. As expected, 
a steeper size distribution reduce the concentration but only 
marginally, so that the dust-to-gas ratio in the interarm re¬ 
gions is slightly enhanced with respect to the previous case. 


3.2 Observational predictions 

In this section, we present results from simulated ALMA 
and HiCIAO observations of the sky models computed via 
RADMC-3D simulations starting from the disc models ob¬ 
tained from hydro simulations. We translate hydro simula¬ 
tions into model observations using the method explained in 
Sect. 2.4. Once having computed the expected full-resolution 
emission maps, we simulate realistic ALMA and HiCIAO im¬ 
ages at (sub-)mm and near-infrared wavelengths. The syn¬ 
thetic images are analysed in terms of the ability to spatially 
resolve the spiral features of our models with the aim to 
test the detectability of gravitationally-induced spiral den¬ 
sity waves. 


3.S.i ALMA simulated observations 

The resulting emission maps at (sub-)mm wavelength 
are used as input sky models for the ALMA simulator (ver¬ 
sion 4.1) in order to produce realistic ALMA observations of 
self-gravitating disc models. The selection of observing pa¬ 
rameters is chosen to ensure enough spatial resolution and 
sensitivity to resolve and detect spiral arms in the inter¬ 
mediate and external region of the discs. Assuming a per¬ 
fect calibration of the visibility measurements, the computa¬ 
tion of ALMA simulated observations include the effects of 
thermal noise from the receivers and atmosphere. The sig¬ 
nal atmospheric corruption effects into the visibilities mea¬ 
surements is evaluated taking into account the atmospheric 
transparency computed using the Atmospheric Transmission 
at Microwaves (ATM) code (Pardo et al. 2001) starting from 
average site conditions and input values for the ground tem¬ 
perature and the Precipitable Water Vapour (PWV) pro¬ 
vided by the user. All the sources are located in the Ophiu- 
cus star-forming region (~ 140 pc) and are observed with a 
transit duration of one hour. 

First, we focus on a representative disc model obtained 
from a given set of model parameters. In Fig. 7 we show 
ALMA simulated observations at different observing wave¬ 
lengths of the disc models obtained starting from an initial 
dust grain size distribution with exponent g = 3.5 and a 
fragmentation velocity equal to 30 m s“^. The observation 
parameters adopted in the images presented here and the 
measured total fluxes are summarized in Table 1. It can 


Frequency 

(GHz) 

PWV 

(mm) 

beam size 
(arcsec x arcsec) 

Total flux 

(Jy) 

rms 

(mjy/beam) 

220 

2.75 

0.063 X 0.055 

0.71 

0.18 

345 

1.26 

0.070 X 0.063 

0.93 

0.49 

460 

1.26 

0.062 X 0.054 

1.12 

0.84 

680 

0.47 

0.062 X 0.057 

1.31 

1.15 


Table 1. Atmospheric conditions, beam size, total flux and rms 
for the simulated observations shown in Fig. 7. 


be noticed that the spiral structure is readily detectable by 
ALMA over a wide range of wavelengths. Note that relative 
intensity between the spiral arms and the external region 
decreases with increasing frequency. This variation is linked 
to the optical regime of the spiral arms at different frequen¬ 
cies. It can be shown that, while the interarm regions are 
optically thin, the arm regions became optically thick for 
frequencies .Si, 345 GHz. Thus, for u Si 345 GHz, the radia¬ 
tion coming from the central protostar is more likely to be 
captured by the arm regions, producing a fainter emission 
from the outer regions. 

In addition, as already noted by Dipierro et al. (2014), 
while the original density structure of the disc is character¬ 
ized by a relatively large number of arms (see Fig. 1), the 
ALMA image show a spiral pattern with two spiral arms 
since the smaller-sized arms have been smeared out by the 
limited resolution of the observations. 

In Fig. 8 is shown the comparison of two ALMA simu¬ 
lated observations at 460 GHz of the disc models obtained 
from simulations with an initial dust grain size distribution 
with power-law exponent g = 3 (left) and g = 4 (right). 
As expected, the slope of the dust grain distribution affects 
the resulting observation. For g = 4, most of the dust mass 
is concentrated in small grains. This produces an optically 
thick inner region that capture all the radiation coming from 
the central protostar. As a result, due to the remarkable de¬ 
crease of the contrast between arm and interarm region, the 
spiral pattern is only barely with ALMA. To illustrate the 
effects of varying the slope of the size grain distribution, in 
Fig. 9 we show the optical depth maps at 460 GHz of the 
disc models obtained from simulations with an initial dust 
grain size distribution with power-law exponent g = 3 (left) 
and g = 4 (right). The optical depth is evaluated locally 
using the following expression: 





da 


Kjj{a) da , 


(14) 


where K^{a) represents the absorption opacity at frequency 
u for the dust species with size a. It can be noted that, 
for grain size distribution with q — 3, the spiral pattern 
region is optically thin. On the contrary, for steeper grain 
size distribution the optical regime is thick for all the disc. 

Observation at (sub-)millimetre wavelengths can probe 
some dust properties in the disc interior. In detail, by mak¬ 
ing some assumptions on the chemical composition, shape 
and size distribution of the dust grains, it is possible to put 
constraints on the level of grain growth. For an optically thin 
disc in the Rayleigh-Jeans limit, the sub-millimetre spectral 
energy distribution (SED) can be approximated as a power- 
law (X K,^ (X where a is the spectral index. For 
grain size distribution of the form adopted in the present 
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Figure 7. ALMA simulated images of the self-gravitating disc model obtained starting from an initial dust grain size distribution with 
exponent q = 3.5 and a fragmentation velocity equal to 30 m s“^. Observing frequencies are reported at the top of each image. Contours 
are 4, 6, 8, 10 x the corresponding rms noise. The white colour in the filled ellipse in the lower left corner indicates the size of the 
half-power contour of the synthesized beam (see table 1 for details). 



Figure 8. Comparison of ALMA simulated images at 460 GHz of disc models obtained from simulations with an initial dust grain size 
distribution with power-law exponent g = 3 (left) and g = 4 (right) and a fragmentation velocity equal to 30 m s~^. Contours are 4, 6, 
8, 10 X the corresponding rms noise: 0.78 mjy/beam (left) and 0.72 mjy/beam (right). The white colour in the filled ellipse in the lower 
left corner indicates the size of the half-power contour of the synthesized beam: 0.062 arcsec x 0.054 arcsec. 
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Figure 9. Comparison of optical depth maps at 460 GHz of disc models obtained from simulations with an initial dust grain size 
distribution with power-law exponent q = S (left) and g = 4 (right). 
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Figure 10. The /3 parameter computed from (left) the model data and from (right) two simulated ALMA observations at 650 fim and 
3 mm of the disc model obtained starting from an initial dust grain size distribution with exponent q = 3 and a fragmentation velocity 
equal to 30 m s“^. The white colour in the filled ellipse in the lower left corner of the right panel indicates the size of the half-power 
contour of the synthesized beam: 0.054 arcsec X 0.052 arcsec. 


work (see Eq. 6), a good approximation for the dust opacity 
is a power law: oc , where the parameter /? is a very 

good indicator of the level of grain growth. While the typi¬ 
cal value for ISM grains is /Siam = 1-7, several observations 
(e.g. Testi et al. 2001, 2003, Natta et al. 2004, Ricci et al. 
2010) have found that /?disc < Asm, which is naturally inter¬ 
preted in terms of grain growth (e.g. Draine 2006). In the 
context of the present work, we calculate the value of the /3 
parameter by computing the ratio of the ALMA simulated 
images at two different frequencies, with the aim to infer the 
different level of dust concentration of large grains in arm 
and interarm regions. The /? parameter can be obtained as: 


In — In 

In — In 1/2 


(15) 


where D, and 7^2 are the fluxes at frequencies and 1 / 2 , 
respectively. In Fig. 10 we show the map of /3 computed 
from the model data (left) and the j3 parameter calculated 
from two simulated ALMA observations at 650 fim and 3 


mm (right). We consider the disc model obtained starting 
from an initial dust grain size distribution with exponent 
q — 3 and a fragmentation velocity equal to 30 m s“^. At 
wavelengths mentioned above we verify that the disc is opti¬ 
cally thin and the dust temperature is high enough to ensure 
that the disc is in the Rayleigh-Jeans limit. Therefore, the 
variations of /3 can be interpreted as a sign of dust trap¬ 
ping regions. From Fig. 10 it can be noticed that the beta 
index map has a shape that resemble the spiral structure 
of the disc model. The spiral arms, where most of the cen¬ 
timetre sized grains are concentrated due to the radial drift, 
are characterized by a lower value of /3. By contrast, the /3 
parameter is relatively higher in the interarm regions, where 
large grains are depleted due to the dust trapping mecha¬ 
nism. Thus, taking into account all the assumptions about 
the dust grain properties, the dust trapping induced by GI 
produces a migration of large dust grains into spiral arms 
region that can be potentially investigated through multi- 
wavelenght observations. 
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3.2.2 Polarized H-band scattered light images 

In this subsection we show synthetic HiCIAO H-band 
polarized intensity image of the disc model obtained start¬ 
ing from an initial dust grain size distribution with exponent 
q = 3.5 and a fragmentation velocity equal to 30 m s“^. As 
previously mentioned, the polarized intensity images essen¬ 
tially trace the gas disc structure since small dust grains 
are less subject to the radial drift. The calculated scattered 
light image convolved with the measured HiCIAO point 
spread function is shown in Fig. 11. The H-band image re¬ 
veals an asymmetry in brightness between the far (top) and 
near (low) disc regions. This is due to forward scattering of 
starlight by large dust particles in a non-face-on disc (incli¬ 
nation of 30°). In addition, the outer ring-like structure is 
very bright since the density of smaller particles increases 
with radius due to the lower value of the fragmentation 
threshold size in outer disc regions (see Fig. 3). Moreover, 
the significative brightness of the outer ring is also caused 
by the fact that scattered light is mostly determined by the 
curvature of the disc surface that, due to the flaring of the 
disc, is very enhanced in the outer regions. However, due to 
the higher brightness of the inner regions and the instrument 
resolution, the spiral structure Is not readily detected. The 
spiral arms and cavities clearly detected in (sub-)millimetre 
ALMA images (see Fig. 7) are smeared out and are thus 
marginally visible in the HiCIAO H-band image. It is worth 
remarking that the shape of the HiCIAO PSF is charac¬ 
terized by large sidelobes which produce a decrease of the 
dynamic range in the convolved images and, consequently, a 
decrease of the detected contrast between arm and interarm 
regions. 

We also consider a new set of simulations with a gas 
disc model two times bigger. In Fig. 12 we present the syn¬ 
thetic HiCIAO H-band polarized intensity image of this spe¬ 
cific disc model with two different inclinations: 10° (left) and 
30° (right). A remarkable improvement of the detection of 
spiral arms can be noticed with comparison to the model 
previously considered (see Fig. 11), due to the fact that the 
HiCIAO PSF is small enough to preserve spiral features in 
scattering images. As expected, the asymmetry in brightness 
between the far (top) and near (low) disc regions due to the 
forward scattering of starlight by large dust particles be¬ 
comes more prominent with increasing angle of inclination. 
In addition, since the separation between arm and interarm 
regions decreases with increasing angle of inclination, the 
detected spiral structure results sharper for face-on discs. 
Moreover, the contrast between arm and interarm region 
appears lower in scattering emission than in thermal emis¬ 
sion. This is caused mostly by the smaller contrast in dust 
density of smaller particles with respect to larger particles 
between arm and interarm regions. All these effect makes it 
more challenging to observe the spiral structures in inclined 
disc through scattering intensity observations. 


4 SUMMARY AND CONCLUSIONS 

In this paper, by combining hydrodynamical and dust 
evolution models, we investigate the influence of gravita¬ 
tional instabilities on the dynamics of dust grains. Using 
representative grains size distributions, we model the dust 



Relative Right Ascension (arcsec) 


Figure 11. Synthetic HiCIAO H-band polarized intensity image 
of the non-face-on disc model (inclination of 30°) obtained starting 
from an initial dust grain size distribution with exponent q = 3.5 
and a fragmentation velocity equal to 30 m s“^. 

grains dynamics under the action of systematic and ran¬ 
dom motion induced by the gas-dust aerodynamical cou¬ 
pling. Once having computed the resulting 2D dust density 
distribution, we perform 3D Monte Carlo radiative transfer 
simulation in order to make predictions of self-gravitating 
discs for future observations. 

The spiral density waves induced by GI produce a re¬ 
markable effect on the dust dynamics leading to the for¬ 
mation of significant dust overdensities in the spiral arms 
regions. The dynamical effect produced by the occurrence 
of pressure bumps in discs closely depends on the size of the 
particles. We found that, assuming that particles are able 
to grow through collisions with velocities up to 10 and 30 
m s“^, the dust density structure of centimetre sized par¬ 
ticles Is affected by the presence of the spiral overdensities. 
While smaller particles closely trace the overall gaseous spi¬ 
ral structure, larger particles experience the highest concen¬ 
tration in density maxima. The influence of the dust trap¬ 
ping mechanism has been tested by computing the dust- 
to-gas ratio. Since most of the mass is concentrated in the 
largest grains, the dust transport process modifies the radial 
and azimuthal distribution of the dust-to-gas ratio leading 
to a creation of a spiral structure. 

The observational predictions of the resulting mod¬ 
els show that the resolution capabilities and sensitivity of 
ALMA are sufficient to spatially resolve the peculiar spi¬ 
ral structure of gravitationally unstable discs with an ac¬ 
ceptable signal-to-noise ratio for non-face-on discs located 
in the Ophiucus (~ 140 pc) star-forming region. Through 
multi-wavelenght (sub-)millimetre observations, we have in¬ 
vestigated whether the density enhancement of the larger 
dust grains into spiral arm regions produces signatures in 
the spectral index map that can be detected using ALMA. 
We find that the migration of larger grains into spiral arms 
region produces clear variations in the spectral index that 
can be interpreted in terms of different level of grain growth 
between arm and interarm regions. 

In addition, our simulations show that the development 
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Figure 12. Synthetic HiCIAO H-band polarized intensity image of a disc model two times bigger than the model of Fig. 11 (i.e. 
Rout = 300 an) with two different inclinations: 10° (left) and 30° (right). 


of gravitational instability can create strong enough surface 
density perturbation that could be detected in near-infrared 
scattered light with HiCIAO. Therefore, the spiral arms ob¬ 
served to date in protoplanetary disc in polarized scattered 
light (Muto et al. 2012, Garufi et al. 2013, Grady et al. 2013) 
might be the result of spiral density perturbations induced 
by the development of gravitational instabilities. There are 
several scenarios that have been proposed to explain the oc¬ 
currence of spirals In protoplanetary discs (Boccaletti et al. 
2013). The more widely accepted one Is based on tidal inter¬ 
action of stellar encounter or embedded protoplanets. How¬ 
ever, Juhasz et al. (2014) have shown that planets by them¬ 
selves cannot induce density perturbation that can be ob¬ 
served in near-infrared scattered light with current telescope 
for sources at the distance comparable to Ophiucus star¬ 
forming region (~ 140 pc). Therefore, the development of 
gravitational instability can be considered a valid theoreti¬ 
cal explanation of the spiral structure observed to date in 
protostellar disc. Moreover, in contrast to the spiral density 
perturbations due to planet-disc interaction (Juhasz et al. 
2014, Dong et al. 2014), the gravitationally-induced density 
waves has a remarkable effect on the dust dynamics and 
produce clear observational signatures at (sub-)millimetre 
wavelengths. 

The simulations presented here contain a number of 
simplifications that need to be addressed. First, since the 
lifetime of the individual spiral features in self-gravitating 
discs is of the order of the approximation of a station¬ 
ary gas density structure for as long as one outer orbital 
period appears to be simplistic (Clarke and Lodato 2009). 
Moreover, in this work we neglect the grain growth and frag¬ 
mentation processing that occur in the disc. This can be 
considered an acceptable approximation since the timescale 
of our simulations is shorter than the typical grain growth 
timescale (Brauer et al. 2008, Okuzumi et al. 2011) in the 
intermediate and outer part of the disc. In addition, it is ex¬ 
pected that the occurrence of long-lived non-axisymmetric 
structures in the disc affects the azimuthal dynamics of dust 
grains leading to an additional systematic drift motion of 


the dust particles towards the pressure maximum (Naka- 
gawa et al. 1986, Birnstiel et al. 2013). Assuming azimuthal 
and radial trapping simultaneously would lead to stronger 
concentrations of dust in pressure traps. The influence of 
time-dependent density perturbations induced by GI and 
the drift caused by azimuthal pressure gradients will be the 
topic of future works. 


ACKNOWLEDGEMENTS 

The authors are grateful to C. P. Dullemond for mak¬ 
ing RADMG-3D available. We thank the referee for use¬ 
ful suggestions. We also thank Attila Juhasz for support 
in solving some problems with RADMC-3D. GD acknowl¬ 
edges Leiden Observatory for hosting and providing CPU 
time in the Sterrewacht workstations. GD and GL acknowl¬ 
edge financial support from PRIN MIUR 2010-2011, project 
“The Ghemical and dynamical evolution of the Milky Way 
and Local Group Galaxies”, prot. 2010LY5N2T. PP is sup¬ 
ported by Koninklijke Nederlandse Akademie van Weten- 
schappen (KNAW) professor prize to Ewine van Dishoeck. 
This work was partly supported by the Italian Ministero 
deiristruzione, Universita e Ricerca through the grant Pro- 
getti Premiali 2012 - iALMA (CUP C52I13000140001). 


REFERENCES 

Birnstiel, T., Dullemond, C. P. and Brauer, F. (2010), A&A 
513, A79. 

Birnstiel, T., Dullemond, C. P. and Plnilla, P. (2013), A&A 

550, L8. 

Blum, J. (2006), Advances in Physics 55, 881. 

Blum, J. and Wurm, G. (2008), ARA&A 46, 21. 

Blum, J., Wurm, G., Kempf, S. et al. (2000), Phys. Rev. 
Lett. 85, 2426. 

Boccaletti, M. et al. (2013), A&A 560, A20. 

Bohren, G. F. and Huffman, D. R. (1983), Absorption and 
Scattering of Light by Small Particles, Wiley. 


© 0000 RAS, MNRAS 000, 000-000 







14 Giovanni Dipierro, Paola Pinilla, Giuseppe Lodato and Leonardo Testi 


Brandt, T. D., McElwain, M. W. et al. (2013), ApJ 

764, 183. 

Brauer, F., Dullemond, C. P. and Henning, T. (2008), A&A 
480, 859. 

Christiaens, V., Casassus, S., Perez, S., van der Plas, G. 
and Menard, F. (2014), ApJ 785, L12. 

Clarke, C. J. and Lodato, G. (2009), MNRAS 398, L6. 
Cossins, P., Lodato, G. and Clarke, C. J. (2009), MNRAS 
393, 1157-1173. 

Cossins, P., Lodato, G. and Testi, L. (2010), MNRAS 
407, 181-188. 

Dipierro, G., Lodato, G., Testi, L. and de Gregorio- 
Monsalvo, 1. (2014), MNRAS 444, 1919-1929. 

Dong, R., Zhu, Z. and Whitney, B. (2014), arXiv:1411.6063 

Draine, B. T. (2006), ApJ 636, 1114-1120. 

Dullemond, C. P., Hollenbach, D., Kamp, 1. and D’Alessio, 
P. (2007), Protostars and Planets V pp. 555-572. 

Durisen, R. H. et al. (2007), Protostar and Planets V 
pp. 607-622. 

Forgan, D., Rice, W. K. M., Cossins, P. and Lodato, G. 
(2011), MNRAS 410, 994. 

Fukagawa, M., Tamura, M. et al. (2006), ApJ 636, L153. 
Gammie, C. F. (2001), ‘Nonlinear outcome of gravitational 
instability in cooling, gaseous disks’, ApJ 553, 174-183. 
Garufi, A., Quanz, S. P., Avenhaus, H. et al. (2013), A&A 
560, A105. 

Gibbons, P. G., Mamatsashvili, G. R. and Rice, W. K. M. 
(2014), MNRAS 442, 361. 

Gibbons, P. G., Rice, W. K. M. and Mamatsashvili, G. R. 
(2012), MNRAS 426, 1444. 

Grady, C. A. M., Muto, T., Hashimoto, J., Fukagawa, M. 
and Currie, T. (2013), ApJ 762, 48. 

Guilloteau, S., Dutrey, A., Pietu, V. and Boehler, Y. 
(2011), A&A 529, A105. 

Gundlach, B. and Blum, J. (2015), ApJ 798, 34. 
Haghighipour, N. and Boss, A. P. (2003a), ApJ 583, 996- 
1003. 

Haghighipour, N. and Boss, A. P. (20036), ApJ 598, 1301- 
1311. 

Johansen, A., Klahr, H. and Henning, T. (2011), A&A 
529, A62. 

Johansen, A., Youdin, A. and Klahr, H. (2009), ApJ 
697, 1269. 

Juhasz, A., Benisty, M. et al. (2014), arXiv:1412.3412 . 
Klahr, H. and Henning, T. (1997), Icarus 128, 213. 

Laibe, G., Gonzales, J. F. and Maddison, S. T. (2012), 
A&A 537, A61. 

Lodato, G. and Rice, W. K. M. (2004), MNRAS 351, 630- 
642. 

Lodato, G. and Rice, W. K. M. (2005), MNRAS358, 1489- 
1500. 

Mathis, J. S., Rumpl, W. and Nordsiecj, K. H. (1977), ApJ 
217, 425-433. 

Miotello, A. et al. (2014), A&A 567, A32. 

Mulders, G. D., Min, M., Dominik, C. et al. (2013), A&A 
549, A112. 

Muto, T., Grady, C. A., Hashimoto, J., Fukagawa, M. and 
Currie, T. (2012), ApJ 748, L22. 

Nakagawa, Y., Sekiya, M. and Hayashi, C. (1986), Icarus 
67, 375. 


Natta, A., Testi, L., Neri, R., Shepherd, D. S. and Wilner, 
D. J. (2004), A&A 416, 179-186. 

Natta, A. et al. (2007), Protostar and Planets V p. 767. 

Okuzumi, S. (2009), ApJ 698, 1122. 

Okuzumi, S., Tanaka, H., Takeuchi, T. and Sakagami, M. 
(2011), ApJ 731, 96. 

Ormel, C. W. and Cuzzi, J. N. (2007), A&A 466, 413. 

Pardo, J. R., Cernicharo, J. and Serabyn, E. (2001), ITAP 
49, 1683. 

Perez, L. M. et al. (2012), ApJ 760, L17. 

Pinilla, P., Birnstiel, T., Ricci, L., Dullemond, C. P., Uribe, 
A. L., Testi, L. and Natta, A. (2012), A&A 538, A114. 

Poppe, T., Blum, J. and Henning, T. (1999), Adv. Space 
Res. 23, 1197. 

Ricci, L., Testi, L., Natta, A., Neri, R., Cabrit, S. and Her- 
czeg, G. J. (2010), A&A 512, A15. 

Rice, W. K. M., Lodato, G. and Armitage, P. J. (2005), 
MNRAS 364, L56-L60. 

Rice, W. K. M., Lodato, G., Pringle, J. E., Armitage, P. J. 
and Bonnell, 1. A. (2004), MNRAS 355, 543-552. 

Rice, W. K. M., Lodato, G., Pringle, J. E., Armitage, P. J. 
and Bonnell, 1. A. (2006), MNRAS 372, L9-L13. 

Shakura, N. 1. and Sunyaev, R. A. (1973), A&A 24, 337. 

Testi, L., Birnstiel, T., Ricci, L. et al. (2014), Protostars 
and Planets VI, astroph/1402.1354 . 

Testi, L., Natta, A., Shepherd, D. S. and Wilner, D. J. 
(2001), ApJ 554, 1087. 

Testi, L., Natta, A., Shepherd, D. S. and Wilner, D. J. 
(2003), A&A 403, 323. 

van der Marel, N. et al. (2013), Science 340, 1199-1202. 

Wada, K., Tanaka, H., Suyama, T., Kimura, H. and Ya¬ 
mamoto, T. (2009), ApJ 702, 1490. 

Weidenschilling, S. J. (1977), MNRAS 180, 57. 

Williams, J. P. and Cieza, L. A. (2011), Annual Review of 
Astronomy and Astrophysics 49, 67-117. 

Youdin, A. N. and Lithwick, Y. (2007), Icarus 192, 588. 

Youdin, A. and Shu, F. H. (2002), ApJ 580, 494. 

Zsom, A. and Dullemond, C. P. (2008), A&A 489, 931. 


© 0000 RAS, MNRAS 000, 000-000 



